home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / sbsp_aux.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  25.2 KB  |  784 lines
  1. /******************************************************************************
  2. * SBsp-Aux.c - Bspline surface auxilary routines.                  *
  3. *******************************************************************************
  4. * Written by Gershon Elber, July. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /* Define some marcos to make some of the routines below look better. They  */
  13. /* calculate the index of the U, V point of the control mesh in Points.        */
  14. #define DERIVED_SRF(U, V)    CAGD_MESH_UV(DerivedSrf, U, V)
  15. #define RAISED_SRF(U, V)    CAGD_MESH_UV(RaisedSrf, U, V)
  16. #define SRF(U, V)        CAGD_MESH_UV(Srf, U, V)
  17.  
  18. /******************************************************************************
  19. * Given a bspline surface - subdivide it into two at given parametric value.  *
  20. * Returns pointer to first surface in a list of two srfs (subdivided ones).   *
  21. * The subdivision is exact result of evaluating the surface int a curve at t  *
  22. * using the recursive algorithm - the left resulting points is left surface,  *
  23. * and the right resulting points is right surface (left is below t).          *
  24. ******************************************************************************/
  25. CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf, CagdRType t,
  26.                             CagdSrfDirType Dir)
  27. {
  28.     CagdBType
  29.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf);
  30.     int i, j, Row, Col, KVLen, Index1, Index2, Mult,
  31.     LULength, RULength, LVLength, RVLength,
  32.     ULength = Srf -> ULength,
  33.     VLength = Srf -> VLength,
  34.     UOrder = Srf -> UOrder,
  35.     VOrder = Srf -> VOrder,
  36.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  37.     CagdRType *RefKV, **Pts, **LPts, **RPts, *LOnePts, *ROnePts, *OnePts;
  38.     CagdSrfStruct *RSrf, *LSrf;
  39.     BspKnotAlphaCoeffType *A;
  40.  
  41.     switch (Dir) {
  42.     case CAGD_CONST_U_DIR:
  43.         RefKV = Srf -> UKnotVector;
  44.         KVLen = UOrder + ULength;
  45.         Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
  46.         Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
  47.         LSrf = BspSrfNew(Index1 + 1, VLength,
  48.                  UOrder, VOrder, Srf -> PType);
  49.         RSrf = BspSrfNew(ULength - Index2 + UOrder, VLength,
  50.                  UOrder, VOrder, Srf -> PType);
  51.         Mult = UOrder - 1 - (Index2 - Index1 - 1);
  52.  
  53.         /* Update the new knot vectors. */
  54.         CAGD_GEN_COPY(LSrf -> UKnotVector,
  55.               Srf -> UKnotVector,
  56.               sizeof(CagdRType) * (Index1 + 1));
  57.         /* Close the knot vector with multiplicity Order: */
  58.         for (j = Index1 + 1; j <= Index1 + UOrder; j++)
  59.         LSrf -> UKnotVector[j] = t;
  60.         CAGD_GEN_COPY(&RSrf -> UKnotVector[UOrder],
  61.               &Srf -> UKnotVector[Index2],
  62.               sizeof(CagdRType) * (ULength + UOrder - Index2));
  63.         /* Make sure knot vector starts with multiplicity Order: */
  64.         for (j = 0; j < UOrder; j++)
  65.         RSrf -> UKnotVector[j] = t;
  66.         
  67.         /* And copy the other direction knot vectors. */
  68.         CAGD_GEN_COPY(LSrf -> VKnotVector,
  69.               Srf -> VKnotVector,
  70.               sizeof(CagdRType) * (VOrder + VLength));
  71.         CAGD_GEN_COPY(RSrf -> VKnotVector,
  72.               Srf -> VKnotVector,
  73.               sizeof(CagdRType) * (VOrder + VLength));
  74.         break;
  75.     case CAGD_CONST_V_DIR:
  76.         RefKV = Srf -> VKnotVector;
  77.         KVLen = VOrder + VLength;
  78.         Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
  79.         Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
  80.         LSrf = BspSrfNew(ULength, Index1 + 1,
  81.                  UOrder, VOrder, Srf -> PType);
  82.         RSrf = BspSrfNew(ULength, VLength - Index2 + VOrder,
  83.                  UOrder, VOrder, Srf -> PType);
  84.         Mult = VOrder - 1 - (Index2 - Index1 - 1);
  85.  
  86.         /* Update the new knot vectors. */
  87.         CAGD_GEN_COPY(LSrf -> VKnotVector,
  88.               Srf -> VKnotVector,
  89.               sizeof(CagdRType) * (Index1 + 1));
  90.         /* Close the knot vector with multiplicity Order: */
  91.         for (j = Index1 + 1; j <= Index1 + VOrder; j++)
  92.         LSrf -> VKnotVector[j] = t;
  93.         CAGD_GEN_COPY(&RSrf -> VKnotVector[VOrder],
  94.               &Srf -> VKnotVector[Index2],
  95.               sizeof(CagdRType) * (VLength + VOrder - Index2));
  96.         /* Make sure knot vector starts with multiplicity Order: */
  97.         for (j = 0; j < VOrder; j++)
  98.         RSrf -> VKnotVector[j] = t;
  99.         
  100.         /* And copy the other direction knot vectors. */
  101.         CAGD_GEN_COPY(LSrf -> UKnotVector,
  102.               Srf -> UKnotVector,
  103.               sizeof(CagdRType) * (UOrder + ULength));
  104.         CAGD_GEN_COPY(RSrf -> UKnotVector,
  105.               Srf -> UKnotVector,
  106.               sizeof(CagdRType) * (UOrder + ULength));
  107.         break;
  108.     default:
  109.         Mult = 1;
  110.         RefKV = NULL;
  111.         LSrf = RSrf = NULL;
  112.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  113.         break;
  114.     }
  115.  
  116.     Pts = Srf -> Points;
  117.     LPts = LSrf -> Points;
  118.     RPts = RSrf -> Points;
  119.     LULength = LSrf -> ULength,
  120.     RULength = RSrf -> ULength;
  121.     LVLength = LSrf -> VLength,
  122.     RVLength = RSrf -> VLength;
  123.  
  124.     switch (Dir) {
  125.     case CAGD_CONST_U_DIR:
  126.         /* Compute the Alpha refinement matrix. */
  127.         if (Mult > 0) {
  128.         CagdRType
  129.             *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  130.  
  131.         for (i = 0; i < Mult; i++)
  132.             NewKV[i] = t;
  133.         A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength,
  134.                                 NewKV, Mult);
  135.         IritFree((VoidPtr) NewKV);
  136.         }
  137.         else
  138.         A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength, NULL, 0);
  139.  
  140.         /* Now work on the two new surfaces meshes. */
  141.  
  142.         /* Note that Mult can be negative in cases where original       */
  143.         /* multiplicity was order or more and we need to compensate     */
  144.         /* here, since Alpha matrix will be just a unit matrix then.    */
  145.         Mult = Mult >= 0 ? 0 : -Mult;
  146.  
  147.         for (Row = 0; Row < VLength; Row++) {
  148.             /* Blend Srf into LSrf. */
  149.             for (j = IsNotRational; j <= MaxCoord; j++) {
  150.             LOnePts = &LPts[j][Row * LULength];
  151.             OnePts = &Pts[j][Row * ULength];
  152.             for (i = 0; i < LULength; i++, LOnePts++)
  153.                 CAGD_ALPHA_BLEND(A, i, OnePts, LOnePts);
  154.         }
  155.  
  156.             /* Blend Srf into LSrf. */
  157.         for (j = IsNotRational; j <= MaxCoord; j++) {
  158.             ROnePts = &RPts[j][Row * RULength];
  159.             OnePts = &Pts[j][Row * ULength];
  160.             for (i = LSrf -> ULength - 1 + Mult;
  161.              i < LSrf -> ULength + RSrf -> ULength - 1 + Mult;
  162.              i++, ROnePts++)
  163.             CAGD_ALPHA_BLEND(A, i, OnePts, ROnePts);
  164.         }
  165.         }
  166.  
  167.         BspKnotFreeAlphaCoef(A);
  168.         break;
  169.     case CAGD_CONST_V_DIR:
  170.         /* Compute the Alpha refinement matrix. */
  171.         if (Mult > 0) {
  172.         CagdRType
  173.             *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  174.  
  175.         for (i = 0; i < Mult; i++)
  176.             NewKV[i] = t;
  177.         A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength,
  178.                           NewKV, Mult);
  179.         IritFree((VoidPtr) NewKV);
  180.         }
  181.         else
  182.         A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength, NULL, 0);
  183.  
  184.         /* Now work on the two new surfaces meshes. */
  185.  
  186.         /* Note that Mult can be negative in cases where original       */
  187.         /* multiplicity was order or more and we need to compensate     */
  188.         /* here, since Alpha matrix will be just a unit matrix then.    */
  189.         Mult = Mult >= 0 ? 0 : -Mult;
  190.  
  191.         for (Col = 0; Col < ULength; Col++) {
  192.         LULength = LSrf -> ULength;
  193.         RULength = RSrf -> ULength;
  194.  
  195.         /* Blend Srf into LSrf. */
  196.         for (j = IsNotRational; j <= MaxCoord; j++) {
  197.             LOnePts = &LPts[j][Col];
  198.             OnePts = &Pts[j][Col];
  199.             for (i = 0;
  200.                  i < LVLength;
  201.                  i++, LOnePts += LULength)
  202.                 CAGD_ALPHA_BLEND_STEP(A, i, OnePts, LOnePts, ULength);
  203.         }
  204.  
  205.         /* Blend Srf into RSrf. */
  206.         for (j = IsNotRational; j <= MaxCoord; j++) {
  207.             ROnePts = &RPts[j][Col];
  208.             OnePts = &Pts[j][Col];
  209.             for (i = LVLength - 1 + Mult;
  210.              i < LVLength + RVLength - 1 + Mult;
  211.              i++, ROnePts += RULength)
  212.             CAGD_ALPHA_BLEND_STEP(A, i, OnePts, ROnePts, ULength);
  213.         }
  214.         }
  215.  
  216.         BspKnotFreeAlphaCoef(A);
  217.         break;
  218.     default:
  219.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  220.         break;
  221.     }
  222.  
  223.     BspKnotMakeRobustKV(RSrf -> UKnotVector,
  224.             RSrf -> UOrder + RSrf -> ULength);
  225.     BspKnotMakeRobustKV(RSrf -> VKnotVector,
  226.             RSrf -> VOrder + RSrf -> VLength);
  227.  
  228.     BspKnotMakeRobustKV(LSrf -> UKnotVector,
  229.             LSrf -> UOrder + LSrf -> ULength);
  230.     BspKnotMakeRobustKV(LSrf -> VKnotVector,
  231.             LSrf -> VOrder + LSrf -> VLength);
  232.  
  233.     LSrf -> Pnext = RSrf;
  234.     return LSrf;
  235. }
  236.  
  237. /******************************************************************************
  238. *  Insert n knot all with the value t in direction Dir. In no case will the   *
  239. * multiplicity of knot be greater or equal to the curve order.              *
  240. ******************************************************************************/
  241. CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *BspSrf, CagdSrfDirType Dir,
  242.                                 CagdRType t, int n)
  243. {
  244.     int i, CrntMult, Mult;
  245.     CagdSrfStruct *RefinedSrf;
  246.  
  247.     switch (Dir) {
  248.     case CAGD_CONST_U_DIR:
  249.         CrntMult = BspKnotFindMult(BspSrf -> UKnotVector, BspSrf -> UOrder,
  250.                              BspSrf -> ULength, t),
  251.         Mult = MIN(n, BspSrf -> UOrder - CrntMult - 1);
  252.         break;
  253.     case CAGD_CONST_V_DIR:
  254.         CrntMult = BspKnotFindMult(BspSrf -> VKnotVector, BspSrf -> VOrder,
  255.                              BspSrf -> VLength, t),
  256.         Mult = MIN(n, BspSrf -> VOrder - CrntMult - 1);
  257.         break;
  258.     default:
  259.         Mult = 0;
  260.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  261.         break;
  262.     }
  263.  
  264.     if (Mult > 0) {
  265.     CagdRType
  266.         *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  267.  
  268.     for (i = 0; i < Mult; i++)
  269.         NewKV[i] = t;
  270.  
  271.     RefinedSrf = BspSrfKnotInsertNDiff(BspSrf, Dir, FALSE, NewKV, Mult);
  272.  
  273.     IritFree((VoidPtr) NewKV);
  274.     }
  275.     else {
  276.     RefinedSrf = CagdSrfCopy(BspSrf);
  277.     }
  278.  
  279.     return RefinedSrf;
  280. }
  281.  
  282. /******************************************************************************
  283. *  Insert n knot with different values as defined by t. If however Replace is *
  284. * TRUE, the knot are simply replacing the current ones.                  *
  285. ******************************************************************************/
  286. CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *Srf, CagdSrfDirType Dir,
  287.                           int Replace, CagdRType *t, int n)
  288. {
  289.     CagdBType
  290.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf);
  291.     int i, Row, Col,
  292.     ULength = Srf -> ULength,
  293.     VLength = Srf -> VLength,
  294.     UOrder = Srf -> UOrder,
  295.     VOrder = Srf -> VOrder,
  296.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  297.     CagdSrfStruct
  298.     *RefSrf = NULL;
  299.  
  300.     if (Replace) {
  301.     for (i = 1; i < n; i++)
  302.         if (t[i] < t[i - 1])
  303.         FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  304.  
  305.         switch (Dir) {
  306.         case CAGD_CONST_U_DIR:
  307.         if (Srf -> UOrder + Srf -> ULength != n)
  308.             FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  309.  
  310.         RefSrf = CagdSrfCopy(Srf);
  311.         for (i = 0; i < n; i++)
  312.             RefSrf -> UKnotVector[i] = *t++;
  313.         break;
  314.         case CAGD_CONST_V_DIR:
  315.         if (Srf -> VOrder + Srf -> VLength != n)
  316.             FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  317.  
  318.         RefSrf = CagdSrfCopy(Srf);
  319.         for (i = 0; i < n; i++)
  320.             RefSrf -> VKnotVector[i] = *t++;
  321.         break;
  322.         default:
  323.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  324.         break;
  325.     }
  326.     }
  327.     else if (n == 0) {
  328.     RefSrf = CagdSrfCopy(Srf);
  329.     }
  330.     else {
  331.     int j, LengthKVt, RULength, RVLength;
  332.     BspKnotAlphaCoeffType *A;
  333.     CagdRType *MergedKVt,
  334.         *UKnotVector = Srf -> UKnotVector,
  335.         *VKnotVector = Srf -> VKnotVector;
  336.  
  337.     switch (Dir) {
  338.         case CAGD_CONST_U_DIR:
  339.         /* Compute the Alpha refinement matrix. */
  340.         MergedKVt = BspKnotMergeTwo(UKnotVector, ULength + UOrder,
  341.                         t, n, 0, &LengthKVt);
  342.         A = BspKnotEvalAlphaCoef(UOrder, UKnotVector, ULength,
  343.                      MergedKVt, LengthKVt - UOrder);
  344.  
  345.             RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder,
  346.                                   Srf -> PType);
  347.         IritFree((VoidPtr) RefSrf -> UKnotVector);
  348.         IritFree((VoidPtr) RefSrf -> VKnotVector);
  349.         RefSrf -> UKnotVector = MergedKVt;
  350.         RefSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
  351.                     Srf -> VLength + Srf -> VOrder);
  352.  
  353.         RULength = RefSrf -> ULength;
  354.  
  355.         /* Update the control mesh */
  356.         for (Row = 0; Row < VLength; Row++) {
  357.             for (j = IsNotRational; j <= MaxCoord; j++) {
  358.             CagdRType
  359.                 *ROnePts = &RefSrf -> Points[j][Row * RULength],
  360.                 *OnePts = &Srf -> Points[j][Row * ULength];
  361.  
  362.             for (i = 0; i < RULength; i++, ROnePts++)
  363.                 CAGD_ALPHA_BLEND( A, i, OnePts, ROnePts );
  364.             }
  365.         }
  366.  
  367.         BspKnotFreeAlphaCoef(A);
  368.         break;
  369.         case CAGD_CONST_V_DIR:
  370.         /* Compute the Alpha refinement matrix. */
  371.         MergedKVt = BspKnotMergeTwo(VKnotVector, VLength + VOrder,
  372.                         t, n, 0, &LengthKVt);
  373.         A = BspKnotEvalAlphaCoef(VOrder, VKnotVector, VLength,
  374.                      MergedKVt, LengthKVt - VOrder);
  375.  
  376.             RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder,
  377.                                   Srf -> PType);
  378.         IritFree((VoidPtr) RefSrf -> UKnotVector);
  379.         IritFree((VoidPtr) RefSrf -> VKnotVector);
  380.         RefSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
  381.                     Srf -> ULength + Srf -> UOrder);
  382.         RefSrf -> VKnotVector = MergedKVt;
  383.  
  384.         RULength = RefSrf -> ULength;
  385.         RVLength = RefSrf -> VLength;
  386.  
  387.         /* Update the control mesh */
  388.         for (Col = 0; Col < ULength; Col++) {
  389.             for (j = IsNotRational; j <= MaxCoord; j++) {
  390.             CagdRType
  391.                 *ROnePts = &RefSrf -> Points[j][Col],
  392.                 *OnePts = &Srf -> Points[j][Col];
  393.  
  394.             for (i = 0; i < RVLength; i++, ROnePts += RULength)
  395.                 CAGD_ALPHA_BLEND_STEP(A, i, OnePts,
  396.                           ROnePts, ULength);
  397.             }
  398.         }
  399.  
  400.         BspKnotFreeAlphaCoef(A);
  401.         break;
  402.         default:
  403.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  404.         break;
  405.     }
  406.     }
  407.  
  408.     BspKnotMakeRobustKV(RefSrf -> UKnotVector,
  409.             RefSrf -> UOrder + RefSrf -> ULength);
  410.     BspKnotMakeRobustKV(RefSrf -> VKnotVector,
  411.             RefSrf -> VOrder + RefSrf -> VLength);
  412.  
  413.     return RefSrf;
  414. }
  415.  
  416. /******************************************************************************
  417. * Return a new surface, identical to the original but with one degree higher  *
  418. * in the given direction.                              *
  419. ******************************************************************************/
  420. CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir)
  421. {
  422.     CagdBType
  423.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  424.     int i, i2, j, RaisedLen, Row, Col,
  425.     Order = Dir == CAGD_CONST_V_DIR ? Srf -> UOrder : Srf -> VOrder,
  426.     Length = Dir == CAGD_CONST_V_DIR ? Srf -> ULength : Srf -> VLength,
  427.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  428.     CagdSrfStruct
  429.     *RaisedSrf = NULL;
  430.  
  431.     if (Order > 2) {
  432.     CagdSrfStruct *UnitSrf;
  433.     int UKvLen1 = Srf -> UOrder + Srf -> ULength - 1,
  434.         VKvLen1 = Srf -> VOrder + Srf -> VLength - 1;
  435.     CagdRType
  436.         *UKv = Srf -> UKnotVector,
  437.         *VKv = Srf -> VKnotVector;
  438.  
  439.     /* Degree raise by multiplying by a constant 1 linear surface in the */
  440.     /* raised direction and constant 1 constant surface in the other.    */
  441.  
  442.     switch (Dir) {
  443.         case CAGD_CONST_U_DIR:
  444.         UnitSrf = BspSrfNew(1, 2, 1, 2,
  445.                     CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  446.         for (i = 0; i < 2; i++)
  447.             UnitSrf -> UKnotVector[i] = i > 0 ? UKv[UKvLen1] : UKv[0];
  448.         for (i = 0; i < 4; i++)
  449.             UnitSrf -> VKnotVector[i] = i > 1 ? VKv[VKvLen1] : VKv[0];
  450.         break;
  451.         case CAGD_CONST_V_DIR:
  452.         UnitSrf = BspSrfNew(2, 1, 2, 1,
  453.                     CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  454.         for (i = 0; i < 4; i++)
  455.             UnitSrf -> UKnotVector[i] = i > 1 ? UKv[UKvLen1] : UKv[0];
  456.         for (i = 0; i < 2; i++)
  457.             UnitSrf -> VKnotVector[i] = i > 0 ? VKv[VKvLen1] : VKv[0];
  458.         break;
  459.         default:
  460.         UnitSrf = NULL;
  461.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  462.         break;
  463.     }
  464.     for (i = 1; i <= MaxCoord; i++)
  465.         UnitSrf -> Points[i][0] = UnitSrf -> Points[i][1] = 1.0;
  466.  
  467.     RaisedSrf = BspSrfMult(Srf, UnitSrf);
  468.  
  469.     CagdSrfFree(UnitSrf);
  470.  
  471.     return RaisedSrf;
  472.     }
  473.  
  474.     /* If surface is linear, degree raising means basically to increase the  */
  475.     /* knot multiplicity of each segment by one and add a middle point for   */
  476.     /* each such segment.                             */
  477.     RaisedLen = Length * 2 - 1;
  478.  
  479.     switch (Dir) {
  480.     case CAGD_CONST_U_DIR:
  481.         RaisedSrf = BspSrfNew(Srf -> ULength, RaisedLen,
  482.                   Srf -> UOrder, Order + 1, Srf -> PType);
  483.  
  484.         /* Update the knot vectors. */
  485.         CAGD_GEN_COPY(RaisedSrf -> UKnotVector,
  486.               Srf -> UKnotVector,
  487.               sizeof(CagdRType) * (Srf -> ULength + Srf -> UOrder));
  488.         for (i = 0; i < 3; i++)
  489.         RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[0];
  490.         for (i = 2, j = 3; i < Length; i++, j += 2)
  491.         RaisedSrf -> VKnotVector[j] = RaisedSrf -> VKnotVector[j + 1] = 
  492.             Srf -> VKnotVector[i];
  493.         for (i = j; i < j + 3; i++)
  494.         RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[Length];
  495.  
  496.         /* Update the mesh. */
  497.                for (Col = 0; Col < Srf -> ULength; Col++) {
  498.         for (j = IsNotRational; j <= MaxCoord; j++)   /* First point. */
  499.             RaisedSrf -> Points[j][RAISED_SRF(Col, 0)] =
  500.             Srf -> Points[j][SRF(Col, 0)];
  501.  
  502.         for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  503.             for (j = IsNotRational; j <= MaxCoord; j++) {
  504.             RaisedSrf -> Points[j][RAISED_SRF(Col, i2)] =
  505.                 Srf -> Points[j][SRF(Col, i - 1)] * 0.5 +
  506.                 Srf -> Points[j][SRF(Col, i)] * 0.5;
  507.             RaisedSrf -> Points[j][RAISED_SRF(Col, i2 + 1)] =
  508.                 Srf -> Points[j][SRF(Col, i)];
  509.             }
  510.         }
  511.         break;
  512.     case CAGD_CONST_V_DIR:
  513.         RaisedSrf = BspSrfNew(RaisedLen, Srf -> VLength,
  514.                   Order + 1, Srf -> VOrder, Srf -> PType);
  515.  
  516.         /* Update the knot vectors. */
  517.         CAGD_GEN_COPY(RaisedSrf -> VKnotVector,
  518.               Srf -> VKnotVector,
  519.               sizeof(CagdRType) * (Srf -> VLength + Srf -> VOrder));
  520.         for (i = 0; i < 3; i++)
  521.         RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[0];
  522.         for (i = 2, j = 3; i < Length; i++, j += 2)
  523.         RaisedSrf -> UKnotVector[j] = RaisedSrf -> UKnotVector[j + 1] = 
  524.             Srf -> UKnotVector[i];
  525.         for (i = j; i < j + 3; i++)
  526.         RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[Length];
  527.  
  528.         /* Update the mesh. */
  529.                for (Row = 0; Row < Srf -> VLength; Row++) {
  530.         for (j = IsNotRational; j <= MaxCoord; j++)   /* First point. */
  531.             RaisedSrf -> Points[j][RAISED_SRF(0, Row)] =
  532.             Srf -> Points[j][SRF(0, Row)];
  533.  
  534.         for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  535.             for (j = IsNotRational; j <= MaxCoord; j++) {
  536.             RaisedSrf -> Points[j][RAISED_SRF(i2, Row)] =
  537.                 Srf -> Points[j][SRF(i - 1, Row)] * 0.5 +
  538.                 Srf -> Points[j][SRF(i, Row)] * 0.5;
  539.             RaisedSrf -> Points[j][RAISED_SRF(i2 + 1, Row)] =
  540.                 Srf -> Points[j][SRF(i, Row)];
  541.             }
  542.         }
  543.         break;
  544.     default:
  545.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  546.         break;
  547.     }
  548.  
  549.     return RaisedSrf;
  550. }
  551.  
  552. /******************************************************************************
  553. * Return a new surface equal to the derived surface in the direction Dir.     *
  554. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
  555. * Q(i) = (k - 1) * P(i+1) - P(i), i = 0 to k-2.                      *
  556. * This is applied to all rows/cols of the surface.                  *
  557. ******************************************************************************/
  558. CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf, CagdSrfDirType Dir)
  559. {
  560.     CagdBType
  561.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  562.     int i, j, Row, Col,
  563.     ULength = Srf -> ULength,
  564.     VLength = Srf -> VLength,
  565.     UOrder = Srf -> UOrder,
  566.     VOrder = Srf -> VOrder,
  567.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  568.     CagdRType **DPoints,
  569.     *UKv = Srf -> UKnotVector,
  570.     *VKv = Srf -> VKnotVector,
  571.     **Points = Srf -> Points;
  572.     CagdSrfStruct
  573.         *DerivedSrf = NULL;
  574.  
  575.     if (!IsNotRational)
  576.     return BspSrfDeriveRational(Srf, Dir);
  577.  
  578.     switch (Dir) {
  579.     case CAGD_CONST_V_DIR:
  580.         if (UOrder < 2)
  581.         FATAL_ERROR(CAGD_ERR_LIN_NO_SUPPORT);
  582.  
  583.         DerivedSrf = BspSrfNew(ULength - 1, VLength,
  584.                    UOrder - 1, VOrder, Srf -> PType);
  585.         CAGD_GEN_COPY(DerivedSrf -> UKnotVector, &UKv[1],
  586.               sizeof(CagdRType) * (ULength + UOrder - 2));
  587.         CAGD_GEN_COPY(DerivedSrf -> VKnotVector, VKv,
  588.               sizeof(CagdRType) * (VLength + VOrder));
  589.         DPoints = DerivedSrf -> Points;
  590.  
  591.                for (Row = 0; Row < VLength; Row++)
  592.         for (i = 0; i < ULength - 1; i++) {
  593.             CagdRType
  594.             Denom = UKv[i + UOrder] - UKv[i + 1];
  595.  
  596.             for (j = IsNotRational; j <= MaxCoord; j++) {
  597.             if (APX_EQ(Denom, 0.0))
  598.                 Denom = INFINITY;
  599.  
  600.             DPoints[j][DERIVED_SRF(i, Row)] = (UOrder - 1) *
  601.                 (Points[j][SRF(i + 1, Row)] -
  602.                  Points[j][SRF(i, Row)]) / Denom;
  603.             }
  604.         }
  605.         break;
  606.     case CAGD_CONST_U_DIR:
  607.         if (VOrder < 2)
  608.         FATAL_ERROR(CAGD_ERR_LIN_NO_SUPPORT);
  609.  
  610.         DerivedSrf = BspSrfNew(ULength, VLength - 1,
  611.                        UOrder, VOrder - 1, Srf -> PType);
  612.         CAGD_GEN_COPY(DerivedSrf -> UKnotVector, UKv,
  613.               sizeof(CagdRType) * (ULength + UOrder));
  614.         CAGD_GEN_COPY(DerivedSrf -> VKnotVector, &VKv[1],
  615.               sizeof(CagdRType) * (VLength + VOrder - 2));
  616.         DPoints = DerivedSrf -> Points;
  617.  
  618.         for (Col = 0; Col < ULength; Col++)
  619.         for (i = 0; i < VLength - 1; i++) {
  620.             CagdRType
  621.             Denom = VKv[i + VOrder] - VKv[i + 1];
  622.  
  623.             for (j = IsNotRational; j <= MaxCoord; j++) {
  624.             if (APX_EQ(Denom, 0.0))
  625.                 Denom = INFINITY;
  626.  
  627.             DPoints[j][DERIVED_SRF(Col, i)] = (VOrder - 1) *
  628.                 (Points[j][SRF(Col, i + 1)] -
  629.                  Points[j][SRF(Col, i)]) / Denom;
  630.             }
  631.         }
  632.         break;
  633.     default:
  634.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  635.         break;
  636.     }
  637.  
  638.     return DerivedSrf;
  639. }
  640.  
  641. /******************************************************************************
  642. * Evaluate the tangent to a surface at a given point and given direction.     *
  643. ******************************************************************************/
  644. CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf, CagdRType u, CagdRType v,
  645.                              CagdSrfDirType Dir)
  646. {
  647.     CagdVecStruct
  648.     *Tangent = NULL;
  649.     CagdCrvStruct *Crv;
  650.  
  651.     switch (Dir) {
  652.     case CAGD_CONST_V_DIR:
  653.         Crv = BspSrfCrvFromSrf(Srf, v, Dir);
  654.         Tangent = BspCrvTangent(Crv, u);
  655.         CagdCrvFree(Crv);
  656.         break;
  657.     case CAGD_CONST_U_DIR:
  658.         Crv = BspSrfCrvFromSrf(Srf, u, Dir);
  659.         Tangent = BspCrvTangent(Crv, v);
  660.         CagdCrvFree(Crv);
  661.         break;
  662.     default:
  663.         FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  664.         break;
  665.     }
  666.  
  667.     return Tangent;
  668. }
  669.  
  670. /******************************************************************************
  671. * Evaluate the normal of a surface at a given point.                  *
  672. * If we fail to compute the normal at given location we try by moving a tad.  *
  673. ******************************************************************************/
  674. CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
  675. {
  676.     static CagdVecStruct Normal;
  677.     CagdVecStruct *V, T1, T2;
  678.     CagdRType UMin, UMax, VMin, VMax;
  679.  
  680.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  681.  
  682.     V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR);
  683.     if (CAGD_LEN_VECTOR(*V) < EPSILON)
  684.     V = BspSrfTangent(Srf,
  685.               u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
  686.               v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
  687.               CAGD_CONST_U_DIR);
  688.     CAGD_COPY_VECTOR(T1, *V);
  689.  
  690.     V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR);
  691.     if (CAGD_LEN_VECTOR(*V) < EPSILON)
  692.     V = BspSrfTangent(Srf,
  693.               u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
  694.               v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
  695.               CAGD_CONST_V_DIR);
  696.     CAGD_COPY_VECTOR(T2, *V);
  697.  
  698.     /* The normal is the cross product of T1 and T2: */
  699.     Normal.Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
  700.     Normal.Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
  701.     Normal.Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];
  702.  
  703.     CAGD_NORMALIZE_VECTOR(Normal);            /* Normalize the vector. */
  704.  
  705.     return &Normal;
  706. }
  707.  
  708. /******************************************************************************
  709. * Evaluate the normals of a surface at a mesh defined by subdividing the      *
  710. * parametric space into a grid of size UFineNess by VFineNess.              *
  711. *   The normals are saved in a linear CagdVecStruct vector which is allocated *
  712. * dynamically. Data is saved u inc. first.                      *
  713. * This routine is much faster than evaluating normal per each point.          *
  714. ******************************************************************************/
  715. CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf, int UFineNess,
  716.                                 int VFineNess)
  717. {
  718.     int i, j;
  719.     CagdRType u, v, UMin, UMax, VMin, VMax;
  720.     CagdVecStruct *Normals, *NPtr, *T, T1, T2;
  721.     CagdCrvStruct **UCurves, **VCurves, *UCrv, *VCrv;
  722.  
  723.     BspSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  724.     Normals = (CagdVecStruct *) IritMalloc(sizeof(CagdVecStruct) * UFineNess *
  725.                                    VFineNess);
  726.  
  727.     UCurves = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) *
  728.                                    UFineNess);
  729.     VCurves = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) *
  730.                                    VFineNess);
  731.  
  732.     UFineNess--;
  733.     VFineNess--;
  734.     for (i = 0; i <= UFineNess; i++)            /* Prepare Iso U curves. */
  735.     UCurves[i] = BspSrfCrvFromSrf(Srf,
  736.                 UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess,
  737.             CAGD_CONST_U_DIR);
  738.     for (j = 0; j <= VFineNess; j++)                /* Prepare Iso V curves. */
  739.     VCurves[j] = BspSrfCrvFromSrf(Srf,
  740.             VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess,
  741.             CAGD_CONST_V_DIR);
  742.  
  743.     NPtr = Normals;
  744.     for (i = 0; i <= UFineNess; i++) {
  745.     UCrv = UCurves[i];
  746.     u = UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess;
  747.  
  748.     for (j = 0; j <= VFineNess; j++) {
  749.         VCrv = VCurves[j];
  750.         v = VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess;
  751.  
  752.         /* We need to copy the tangents as BspCrvTangent save it in as   */
  753.         /* static so the second call will overwrite first call value.    */
  754.         /* If we fail to compute the tangent, we try adjacent isoline.   */
  755.         T = BspCrvTangent(UCrv, v);
  756.         if (CAGD_LEN_VECTOR(*T) < EPSILON)
  757.         T = BspCrvTangent(UCurves[i == 0 ? i + 1 : i - 1], v);
  758.         CAGD_COPY_VECTOR(T1, *T);
  759.  
  760.         T = BspCrvTangent(VCrv, u);
  761.         if (CAGD_LEN_VECTOR(*T) < EPSILON)
  762.         T = BspCrvTangent(VCurves[j == 0 ? j + 1 : j - 1], u);
  763.         CAGD_COPY_VECTOR(T2, *T);
  764.  
  765.         /* The normal is the cross product of T1 and T2: */
  766.         NPtr -> Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
  767.         NPtr -> Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
  768.         NPtr -> Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];
  769.  
  770.         CAGD_NORMALIZE_VECTOR(*NPtr);           /* Normalize the vector. */
  771.         NPtr++;
  772.     }
  773.     }
  774.  
  775.     for (i = 0; i <= UFineNess; i++)
  776.     CagdCrvFree(UCurves[i]);
  777.     IritFree(UCurves);
  778.     for (j = 0; j <= VFineNess; j++)
  779.     CagdCrvFree(VCurves[j]);
  780.     IritFree(VCurves);
  781.  
  782.     return Normals;
  783. }
  784.